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i , Abstract 

We present molecular dynamics simulations of the thermodynamic 
melting transition of a bcc metal, vanadium using the Finnis-Sinclair 
potential. We studied the structural, transport and energetic proper- 
ties of slabs made of 27 atomic layers with a free surface. We investi- 
gated premelting phenomena at the low-index surfaces of vanadium; 
^ . V(lll), V(001), and V(011), finding that as the temperature increases, 

O the V(lll) surface disorders first, then the V(100) surface, while the 

in V(110) surface remains stable up to the melting temperature. Also, as 

. the temperature increases, the disorder spreads from the surface layer 

into the bulk, establishing a thin quasiliquid film in the surface region. 
We conclude that the hierarchy of premelting phenomena is inversely 
proportional to the surface atomic density, being most pronounced for 
the V(lll) surface which has the lowest surface density. 
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Theories of melting [UEHHIZI can be separated into two classes. The first one 
describes the mechanical melting of a homogeneous solid resulting from lat- 
tice instability 0E1 IE] and/or the spontaneous generation of thermal defects 
pi lTul[T^lll3j . The second class treats the thermodynamic melting of hetero- 
geneous solids, which begins at extrinsic defects such as a free surface or an 
internal interface (grain boundaries, voids, etc) [U El HU E3 EH EH UH EE] • 
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Both types of melting have been studied extensively for the case of fee metals. 
In order to determine whether the mechanism of homogeneous melting differs 
if the lattice symmetry is changed, in our first paper [I"""J , we extended these 
studies to a bec metal, vanadium. We used molecular dynamics simulations 
to study the mechanical melting transition in a bulk system, and found that 
the melting transition occurs uniformly throughout the solid at the melting 
point. The shear elastic instability leading to this melting of vanadium oc- 
curred once the molar volume reached a critical value, whether reached by 
heating the solid or by adding defects at a constant temperature [Hj. The 
temperature at which the mechanical melting transition takes place in our 
model is T s = 2500 ±10 K, which is above the melting temperature measured 
experimentally, T m = 2183 K. Thus, the mechanical melting transition can 
be observed only if melting at surfaces is artificially suppressed, since in the 
laboratory it is preempted by the thermodynamic melting transition which 
begins at a free surface. 

Theoretical aspects of melting beginning at the surface have been stud- 
ied by phenomenological theories [T""J CHI 033 1213 I2J] , lattice-gas models [T7] . 
and density functional theory (221- Microscopic descriptions of static and 
dynamic properties of melting phenomena beginning at a surface emerged 
from computer simulations which employ many-body interaction potentials 
derived from the effective-medium j231 I2H EH| and embedded-atom the- 
ories [23 |2E1 I2B EU I2H]- as well as pairwise interatomic interactions in 
the form of Lennard- Jones potentials [3*2*1 I3*3*j . Melting at the surface is 
anisotropic, meaning that certain surfaces of fee metals (Pb(110), Al(110)) 
exhibit premelting,|16. while this phenomenon is not observed for the close- 
packed surfaces (Pb(lll), Al(lll)). Theoretical aspects of this anisotropy 
were discussed by Trayanov and Tosatti [T7j . 

This paper is organized as follows: details of the molecular-dynamics 
simulations are presented in Sec. II. Sec. Ill contains most of the results 
of this study, including the thermodynamic melting temperature and several 
structural, energetic, and transport properties of the surface. Finally, in Sec. 
IV, we summarize our results. 

2 Simulation details 

We model the melting of vanadium with a free surface using molecular dy- 
namics (MD) simulations [3*3 EDj in a canonical ensemble. The many-body 
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Figure 1: Geometry of the sample (V(001) surface). 

interaction potential used in this study was developed by Finnis and Sinclair 
[3~%] (FS). The computational aspects of the algorithm are described in detail 
in our previous article |14j . 

Ideally, a crystal with a surface is a semi-infinite system. We modeled this 
system as a thick slab. This slab is shown in Fig. ^ The bottom 3 atomic 
layers have the positions of the atoms fixed in order to mimic the effect of 
the presence of the infinite bulk of the crystal. Due to the relatively short 
range of the repulsive part of the modified FS potential, 3 fixed layers are 
sufficient to represent the static substrate. On top of those 3 layers there are 
24 layers in which the atoms are free to move. Periodic boundary conditions 
are used along the in-plane (x and y) directions, while the top boundary of 
the slab, along the z direction, is free. In Figure ^ we show a sample where 
the temperature is high enough so that the layers near the free surface are 
already disordered. 

Three different samples with various low-index surfaces were constructed: 
V(001), V(011) and V(lll), whose surface layers are shown in Fig. |21 The 
V(001) and V(lll) samples contain 2700 atoms, initially arranged as a per- 
fect bcc crystal of 10x10x27 unit cells (100 atoms in a layer). The V(011) 
sample contains 2646 atoms arranged as 7x7x27 unit cells (98 atoms in a 
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Figure 2: Geometry of the low-index faces of vanadium (top view). The 
drawings are to scale. 



Each simulation starts from a low-temperature solid. The lattice constant 
of the frozen layers is assigned a value appropriate to the bulk^l- As the 
temperature is raised, this value is adjusted to reflect thermal expansion as 
obtained from bulk simulations |14j. We found that ~ 10 5 integration time 
steps were sufficient in order to reach thermal equilibrium after a temper- 
ature change (a time step dt = 1.05 x 1CT 15 sec). An equilibrium state 
is considered to be achieved when there is no significant temporal variation 
(beyond the statistical fluctuations) of the total energy, layer occupation 
number, structure order parameters, self diffusion coefficients, and a uniform 
profile of kinetic temperature across the sample is observed. Thereafter the 
various structural and transport properties of the system are calculated, ac- 
cumulated and averaged over a long period of 10 7 MD steps. Throughout this 
study, interactive visualization (the AViz program jJT]) was implemented, to 
observe sample disorder and layer mixing (see FigQ). A discussion of this 
visualization can be found in Ref. [4*5] . 

3 Premelting effects on surfaces 

3.1 Thermodynamic melting temperature 

In order to investigate the phenomenon of surface disorder and premelting, we 
first determined the thermodynamic melting point T m of vanadium described 
by the modified FS potential. T m was obtained by means of the method 
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proposed by Lutsko et al. [12] • In this way, we found T m = 2220 ± 10K 
for the V(lll) sample, which has the lowest surface packing density. This 
value is close to the experimental value [13] of T m = 2183/C. For the other 
surfaces, the value of T m comes out slightly higher, but within the uncertainty 
limits of the simulation. The difference could be due to superheating effects 
on the solid-liquid interface just above the thermodynamical melting point 
T m . The largest difference is found for the close packed V(011) surface, with 
T m m 2240 ± 10 K. This dependence of the thermodynamic melting point on 
the crystalline plane forming the surface was also observed in MD simulations 
of fee metals (23 123 123 123 EES] • 

3.2 Surface thermal expansion and amplitude of atomic 
vibration 

The interlayer relaxation and surface thermal expansion was calculated using 
the difference between the average z coordinate of the zth and i + 1th layers: 



where Zj is the ^-coordinate of the atom i, rii is the instantaneous number 
of atoms in the layer i, and the angular brackets denote a time average. 
The distances between the neighboring layers could be determined up to the 
temperatures where surface premelting effects blur the distinction between 
individual layers, although some structure remains visible in the local density 
profiles (see below). 

At low temperatures, the first layer exhibits an inward relaxation, i.e. 
A 1,2 < 0, where we define 



where d i)i+ i is the distance in the ^-direction between the i and i + 1 surface 
layers, and d huik is the distance between two neighboring layers in the bulk. 
In Figures 01 and |U A^+i is plotted versus temperature for the two of the 
low-index faces. 

The effect of the different geometry of the various surfaces is reflected in 
the thermal expansion, i.e. a surface with a lower atomic density expands 
more than a surface which is close packed. As shown in Fig. |3 the thermal 
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Figure 3: Plot of A^+i for the first (squares), second (triangles) and third 
(diamonds) surface layers as a function of temperature for the V(011) surface. 
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Figure 4: Plot of Aj j +1 for the first (squares) and second (triangles) surface 
layers as a function of temperature for the V(lll) surface. 
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expansion of the V(lll) surface layers is significantly larger than the thermal 
expansion of other surfaces and fixed layers (bulk) . In Fig. we plot thermal 
expansion of surface layers (for the V(001) and V(011) slabs) and fixed layers 
on an enlarged scale. The observed "anomalous" thermal expansion of these 
surface layers is a direct consequence of enhanced vibration at the surface, 
where atoms probe the more anharmonic region of the interatomic potential. 
Anomalous thermal expansion was first observed in experiments on Pb(llO) 
by Frenken |15j . 
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Figure 5: Normalized interlayer distances between the first and second layers, 
d 1>2 (T)/d lj2 (T ) at T = WOOK for the 1/(111) (squares), 1/(001) (diamonds), 
1/(011) (triangles) surfaces and the bulk (crosses). Note the anomalous ther- 
mal expansion of the 1/(111) surface. 

Atomic vibration properties at the different faces of vanadium are impor- 
tant for illustrating the onset of anharmonicity at high temperatures. We 
calculate the mean square vibration amplitude (MSVA), < u 2 a >, around the 
equilibrium position in the first surface layer as: 

1 ni /r i2\ 
< U l >= — E \ rUt)- < fi,a(t) >\ ) (3) 
fi l i=l N ' 

where a = x,y, z denotes the spatial direction of motion, rii is the instan- 
taneous number of atoms in the first surface layer, and the angular brackets 
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Figure 6: Normalized interlayer distances between the first and second lay- 
ers, d 1>2 (T)/d 1)2 (To) at T = 1000.R: for the V(001) (diamonds) and 1/(011) 
(triangles) and the bulk (circles). The solid lines are drawn to guide the eye. 

represent an average over time. 

In a harmonic system, < u 2 a > increases linearly with temperature, and 
therefore deviations from linearity are a measure of anharmonicity. The 
mean square amplitudes of vibration shown in Fig. [7|for the V(lll) surface 
demonstrate substantial anharmonicity at elevated temperatures. Figure 
also shows that there is a noticeable difference between the in-plane MSVA, 
< u 2 >, and out-of-plane one, < u 2 z > (we found that < u 2 x >~< u 2 >). 

In-plane MSVA's are larger than out-of-plane amplitudes also for V(001) 
and V(011) surfaces. One reason for this anisotropy could be the different 
nearest neighbour distance along the x and z directions caused by the inward 
relaxation of the surface layer, leading to a net restoring force that is higher 
in the direction perpendicular to the surface. The in-plane components of 
the mean square amplitudes of vibration for the various faces of vanadium 
are shown in Fig. |HJ As may be seen from the figure the in-plane MSVA 
is largest for the V(lll), and smallest for the V(011)surface. Atoms which 
belong to the loose-packed surface V(lll) are less tightly bound than atoms 
of the close-packed V(011) face, and therefore vibrate with larger amplitudes. 
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Figure 7: Mean-square amplitudes of vibration < u 2 z > (filled circles) and 
< u 2 x > (squares) for the first layer of the V(lll) surface as a function of 
temperature. The solid lines are a fit to the data. 
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Figure 8: The calculated in-plane compoments of mean square amplitudes 
of vibration < v? > for the first surface layer of V(lll) (filled circles), V(001) 
(croses) and V(011) (squares) as a function of temperature. The solid lines 
are a fit to the data. 
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3.3 Layer density profiles 



In order to display the structure of the sample along the z direction, perpen- 
dicular to the surface, it is convenient to use the density p(z), defined so that 
p(z)dz is the number of atoms in a slice of thickness dz at z. We represent 
p(z) by a continuous function defined according to Chen et al. [271 12E] 

where Z{ is the z coordinate of atom i, with z = set at the bottom of 
the first layer which is not fixed, and the angular brackets indicate a time 
average. We use Az = 0.1ao/2v3, where ao is the bulk lattice parameter 
at a given temperature. For a system in equilibrium, p(z) can be obtained 
by an average over many different configurations. Plots of p(z) at different 
temperatures are given in Fig. |U]for the V(lll) slab, in Fig.^^for the V(011) 
slab and in Fig. EH for the V(001) slab. 
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Figure 9: Density profile across the V(lll) slab along the z direction at 
various temperatures. 

Premelting effects at the crystal surface can be examined by monitoring 
the layer- by-layer modulation of the density profile of the system, p(z), at 
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Figure 10: Density profile across the V(011) slab along the z direction at 
various temperatures. 
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Figure 11: Density profile across the V(001) slab along the z direction at 
various temperatures. 
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various temperatures up to the melting point T m . At low temperatures p(z) 
consists of a series of well resolved peaks. The atoms are packed in the 
layers with constant density in each layer and virtually no atoms in between 
these layers. As the temperature increases the effective width of each layer 
becomes broader due to the enhanced atomic vibration. The position of 
the density peaks moves to larger values of z due to the thermal expansion. 
At higher temperatures the atomic vibration becomes so large, especially in 
the topmost layers, that the minima of p(z) in between two layers rise to 
non-zero values. In addition, disorder sets in, with atomic migration taking 
place between the layers. Comparing Figs. IMTTl one can see that the V(lll) 
surface crosses over to a state of premelting first, than V(001), while on the 
V(011) surface premelting effects are very small. 

Above some temperature (T pm ~ 1000 K for the V(lll) surface, T pm ~ 1800 
for the V(001) surface, T pm ~ 2200 K for the V(011) surface) the density of 
the topmost layer becomes slightly lower than that of the underlying layers. 
This loss of density is compensated by the appearance of an adlayer. The 
adlayer appears first on the least packed surface V(lll), and then at a higher 
temperature on the V(001) surface. The close packed face V(011) develops 
an adlayer at temperatures very close to T m . This hierarchy indicates that 
the formation energy of structural defects (vacancy-adatom pairs) is differ- 
ent on these surfaces. That hierarchy was also observed in the investigation 
of the surface premelting of fee metals (A1|3D E01, Ni jZ3 |2Hl HE 123 , and 
Cu [2ZII2S1)) where it was found that an adlayer appears first on the least 
packed (011) surface, then at higher temperature on the (001) face, and fi- 
nally on the close packed (111) surface at a temperature close to T m . As the 
temperature increases towards T m the distinction between the surface layers 
becomes blurred, which is consistent with the topmost layers converting into 
a quasiliquid. 

3.4 Layer structure factors 

The preceding section dealt with the onset of disorder along the z direction, 
perpendicular to the surface. In order to monitor the onset of in-plane dis- 
order in each atomic layer, it is useful to define a structure factor (order 
parameter), rji >a , through a Fourier transform of the atomic density which is 
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Figure 12: In plane structure factor, r] x , of the V(lll) slab along the x 
direction vs layer number at several temperatures: T = 40LW (circles), T = 
861W (squares), T = 1100X (diamonds), T = 150LW (triangles down), 
T = 1800AT (triangles up), T = 210LW (triangles left), T = 220LW (triangles 
right), and T = 223LW (crosses). 
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Figure 13: In-plane structure factor, r) x , of the V(011) slab calculated along 
x direction vs layer number at several temperatures: T = lOOO-ft' (circles), 
T = 1200fT (squares), T = 150LW (diamonds), T = 1700K (triangles down), 
T = 1900X (triangles up), T = 2100K (triangles left), T = 2231W (triangles 
right), and T = 2250K (crosses). 
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calculated separately for each layer 



(5) 
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where a = x,y and g a = 2ir/a a e a is a set of reciprocal lattice vectors. a a is 
the distance between nearest neighbors, and rti is the instantaneous number 
of atoms in the layer I. The sum extends over the particles in the layer I, 
and the angular brackets denote averaging over time. 

For an ordered crystalline surface the structure factor is a unity at zero 
temperature. Enhanced vibration and formation of point defects lead to a de- 
crease of the structure factor with increasing temperature. This is illustrated 
in Figs. H21 and Uni where the structure factor along the x direction is plotted 
vs. layer number for the V(lll) and the V(011) slab, respectively. It is seen 
that these effects are most pronounced in the surface region. Vacancies do 
not directly affect the order parameter, since a normalization procedure is 
employed during each measurement by using the instantaneous layer occu- 
pation, rii, of a layer. Nevertheless, vacancies have an indirect effect on the 
order parameter by introducing a localized lattice distortion. 

As Figs. H21 and [T3] show, for the V(lll) slab we found a continuous 
decrease of the in-plane order parameter. In contrast, the close packed V(011) 
sample exhibits a relatively abrupt decrease of the structure factor within less 
than ~ 20K of the melting temperature. 

3.5 Radial distribution function 

The formation of a quasiliquid film can be analyzed by using a plane radial 
distribution function defined as 



where r^y is the component of the r*j — fj parallel to the surface plane, rii 
is the instantaneous number of atoms in layer I, the sum extends over all 
particles in layer I, and the angular brackets denote averaging over time. 

The two-dimensional radial distribution function, p(rii), for several lay- 
ers of the V(lll) and V(011) slabs very close to T m is shown in Figs. El 
and ED As seen from the figures the intra-layer structure in these layers 
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changes gradually from crystalline to liquid-like as one approaches the sur- 
face. Particularly noticeable is the disappearance of the crystalline features 
in p(r\\) that correspond to the second, third and other nearest neighbors. In 
addition to the heights of the peaks, the area under the p(rn) curve changes, 
which reflects the change in the density across the solid-liquid interface. 
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Figure 14: Two-dimensional radial distribution function p(r\\) of the top 
surface layers of V(lll) at temperature T=2230K. The layer n=0 corresponds 
to the adlayer, n=l to the first layer, n=2 to the second one, etc. 

We note that within the adlayer, the probability of finding particles with 
separation beyond the first-neighbor shell is relatively small, indicating a 
tendency for clustering which persists, though to a smaller extent, even to 
T ~ T m . The plane pair-correlation functions for the first layer of the V(001) 
at elevated temperatures are shown in Fig. It is clear that the crystalline 
structure vanishes gradually and the layer becomes quasiliquid as the melting 
point T m is approached. 

We conclude from the analysis of the plane radial distribution functions 
and density profiles that surface premelting begins first on the least packed 
surface V(lll), while changes only appear on the closed packed V(011) sur- 
face at temperatures which are very close to the melting point. 
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Figure 15: Two-dimensional radial distribution function p(r\\) of the surface 
layers of V(011) at temperature T=2230K. The layer n=0 corresponds to the 
adlayer, n=l to the first layer, n=2 to the second one, etc. 
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Figure 16: Two-dimensional radial distribution function p(r\\) of the first 
layer of V(001) at various temperatures. 
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3.6 Layer occupation and point defects 

As seen in the density profiles (Figs. EJ-refdensOOl), the formation of an ad- 
layer on the least packed surface V(lll) begins at around T ~ 800 K. 
At these temperatures, the appearance of an adlayer involves generation 
of vacancies only in the first surface layer, while at higher temperatures 
(T > 1Q00K) vacancies in the underlying layers (the second and third lay- 
ers) begin to appear via promotion of atoms to vacant sites in the first surface 
layer. Atom migration from the deeper layers into the surface layers increases 
significantly as the temperature approaches T m . In contrast to the V(lll) 
sample, adlayer formation and generation of adatom-vacancy pairs at the the 
V(001) surface becomes observable only above T ~ 1400i^ (in practice, all 
adatoms come from the first layer). 

The same effect, namely the appearance of an adlayer at different surfaces 
at successively higher temperatures was observed in computer simulations of 
surface premelting of fee metals [2312312111231211, where the lowest density 

(110) surface of fee metals begins to disorder first, while the close packed 

(111) surface preserves its ordered crystalline structure almost up to T m . 
Knowledge of the equilibrium averaged number of atoms in the adlayer 

allows us to estimate E s , the formation energy of a surface defect (adatom- 
vacancy pair) according to the relation: 

n/N(T = 0) = exp(-E./k B T) (7) 

where n is the adlayer occupation number at a temperature T, and N(T = 0) 
is the number of atoms in a layer at zero temperature. Eq. (7) holds as long 
as the surface structure is preserved, and interactions between defects can 
be neglected. Hence, we used the adlayer occupation data in a temperature 
region in which the concentration of adatom-vacancy pairs is small. The 
natural logarithm of adlayer occupation vs. 1 / k{T is shown in Fig. El for 
the V(lll), V(001) and V(011) surfaces, respectively. The calculated surface 
defect formation energies of the various low-index faces of vanadium are 
tabulated in Table [TJ 

An alternative method of estimating values of E s is to use the contribution 
of defects to anharmonicity, as reflected in the temperature dependence of 
the mean square vibration amplitude < u 2 >. Because the lattice in the 
vicinity of a defect is distorted, the value of < u 2 > will be larger for atoms 
near a defect. The layer averaged contribution to < u 2 > will therefore be 
proportional to the number of defects n, and hence to exp(—E s /kBT). To 
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Figure 17: The natural logarithm of the equilibrium adatom concentration, 
n, as a function of l/k b T on the V(001), V(011) and V(lll) faces (denoted 
by diamonds, triangles and squares, respectively.) The dashed lines are linear 
fits to the data. 

extract E s , we have fitted our < u 2 > data to the form 

< u 2 >= aT + b{T) exp(-E s /k B T) (8) 

where a is a constant, and b(T) is a low order polynomial in T. We found 
that taking b(T) to be first order in T, namely b(T) = bT, was sufficient 
to obtain a good fit. Typical fits were shown as solid lines in Figs. [7| and 
IH1 Values of E s deduced from these fits for the various surfaces are given in 
Table [T] There is good agreement between the values of E s obtained by the 
two methods. 

To summarize this subsection, the formation energy of surface defects 
is lowest at the V(lll) surface, and increases at the V(001) and V(011) 
surfaces. This hierarchy is consistent with the results obtained for various 
faces of copper (fee lattice) by Hakkinen et al. [23] In both cases, the surface 
defect formation energy is largest for the close packed surfaces (V(011) in 
case of a bec lattice, and Cu(lll) in case of a fee lattice), and lowest for the 
least packed surfaces, V(lll) and Cu(Oll), respectively. 
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Table 1: Formation energy, E s , of adatom- vacancy pairs. Top line-Eq.(7), 
second line-Eq.(8). The data for Cu is from Ref. (22] Units are in eV. 



Surface 


(111) 


(001) 


(011) 


Va (bee) 


0.6 ±0.2 eV 


0.67 ±0.03 eV 


2.68 ±0.2 eV 


Va (bee) 


0.53 ±0.05 eV 


0.73 ±0.05 eV 


2.5 ±0.3 eV 


Cu (fee) E s 


1.92 eV 


0.86 eV 


0.39 eV 



3.7 Diffusion coefficients 



Mass transport at surfaces can be characterized using the self diffusion co- 
efficients. These coefficients are found from the particle trajectories, r^(t), 



by calculating the average mean square displacement Rf 



Rl = {-j:inAt + T)-nAT)} 2 ) (9) 

\ ' iei I 

where /i = x, y, z is a coordinate index, the sum includes atoms in the layer 
and the angular brackets denote averaging over time from the origin (r). 
The diffusion coefficients Di ^ are calculated separately for each layer in 
the x, y and z directions, according to the Einstein relation 

R 2 

D lu = lim (10) 

Values of the diffusion coefficients versus layer number are shown in Fig. 
As expected, the mobility of atoms is larger the closer one comes to the sur- 
face. At high enough temperatures, where the surface layers become quasiliq- 
uid, the diffusion coefficient is the same (within the error bars) for all the 
atoms in the region of the surface. These observations correlate with the 
structural variations in the surface region exhibited in the pair correlation 
functions, the structure order parameters, and local density profiles. 

The natural logarithm of the in-plane diffusion coefficients, D\\ = \{D X + 
D y ) , for the first surface layer of the various faces vs. temperature is shown in 
Fig. Ell Within our accuracy, the diffusion can be characterized by a simple 
exponential dependence throughout the whole temperature range, D(T) oc 
exp (—Ed/ksT), with diffusion activation energy, Ed, being independent of 
temperature. Values of E^ for the various low-index faces are given in TableEl 
The diffusion coefficient of the least packed V(lll) face is the largest. The 
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Figure 18: Diffusion coefficients in layers, vs temperature, for the V(lll) 
system along the x-direction ([110]). Squares correspond to the adlayer, 
1 = 0, triangles to the first layer, 1 = 1, diamonds to the second layer, I = 2 
and crosses denote coefficients of diffusion in layer 1 = 3. 

Table 2: Calculated diffusion activation energy and migration energy E m 



Surface 


V(lll) 


V(001) 


V(011) 


E d (eV) 


0.96 ±0.02 eV 


1.56 ±0.03 eV 


3.29 ±0.2 eV 


£ m (eV) 


0.43 ±0.07 eV 


0.87 ±0.08 eV 


0.79 ±0.5 eV 



diffusion coefficient of V(001) is smaller, but still larger than the diffusion 
coefficient of the close packed V(011) surface. In the temperature range where 
the surface region is still solid-like, diffusion takes place mainly by atoms 
exchanging places with vacancies. Writing as a sum of the formation 
energy, E s , and migration energy, E m , and using the values of E s from the 
preceding section, we can obtain values of E m for the various surfaces. These 
are listed in Table El It is seen that the migration energy for the V(lll) 
surface is again the lowest, while for the other two surfaces E m is the same 
within the error bars of our calculation. 
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Figure 19: Plot of natural logarithm of the in-plane diffusion coefficient 
log(D||), for the various faces of vanadium as a function of inverse tempera- 
ture. (D\\) is in units of A 2 /psec. 

4 Conclusions 

We investigated premelting effects at the various surfaces of vanadium using 
molecular dynamics simulations. We determined the thermodynamic melting 
temperature of vanadium described by the FS potential as 2220 ± 10 K 
(V(lll) surface), and studied structural, transport (diffusion) and energetic 
properties (formation energy of surface defects). 

We found that the surface region of 1/(111) begins to disorder first, via 
the generation of vacancy-adatom pairs and a formation of an adlayer at 
temperatures above 1000 K. At higher temperatures, the surface region be- 
comes quasiliquid. This process begins above 1600 K for the V(001) surface. 
At the closest packed V(011) surface, this effect is observed only in close 
proximity of T m . 

The results of our simulations of surface premelting of the bcc metal, 
vanadium, are similar to the results obtained for various fee metals, in the 
sense that the onset of disorder is seen first at the surface with the lowest 
density. 
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